package edu.gmu.atelier;


/**
 * Performs Eigen factorization (a.k.a. decomposition) on the passed matrix.
 * The method involves recursion (i.e. we do not try to determine the 
 * roots of the characteristic equation). This implementation uses 
 * the QRFactor iteratively, QRFactor(A0), A1 = R0 Q0, QRFactor(A01), etc.
 * Taken from Olver and Shakiban Chapter 10.6
 * 
 * @author James H. Pope
 */
public class EigenQRFactor implements EigenFactor
{
    public static final int DEFAULT_MAX_ITERATIONS = 15;
    
    private Matrix a  = null;
    
    private Matrix d  = null; // upper case lambda, eigen values matrix
    
    private Matrix s  = null;
    
    private Matrix q  = null;
    
    private int maxIterations = -1;
    
    private int iterations = -1;
    
    public EigenQRFactor( Matrix a )
    {
        // Default is GRAM_SCHMIDT, may eventually change to HOUSEHOLDER
        this(a, QRFactorFactory.GRAM_SCHMIDT, DEFAULT_MAX_ITERATIONS);
    }
    
    public EigenQRFactor( Matrix a, QRFactorFactory f )
    {
        // Default is GRAM_SCHMIDT, may eventually change to HOUSEHOLDER
        this(a, f, DEFAULT_MAX_ITERATIONS);
    }
    
    public EigenQRFactor( Matrix a, QRFactorFactory f, int maxIterations )
    {
        this.a = a;
        this.maxIterations = maxIterations;
        
        // How many times do we iterate? Currently unknown
        
        /*
         * Peculiar step, note RQ instead of QR, so different A_k.
         * For Symmetric matrices each iteration slowly turns:
         *   R into D with eignen values on the diagonal and zeroes elsewhere
         *   S into the orthonormal eigen vectors
         *   Q into I the identity matrix
         * The end result is intuitive because the "system" settles down
         * and AS = e0q0+e1q1+e2q2 = DS.
         */
        Epsilon e = Epsilon.makeDecimal(15);
        //Matrix i = a.identity();
        // First iteration 0
        QRFactor f0 = f.makeQRFactor(a);
        Matrix qk = f0.getQ();
        Matrix rk = f0.getR();
        Matrix sk = qk;
        Matrix ak_1 = rk.mult(qk);
        // We have already performed one iteration
        int k = 1;
        for( ; k < maxIterations; k++ )
        {
            /*
             * We want to speed up convergence, below method not working?
             * From Strang Page 488, shift before, shift back after
             * Shift by suspected eigenvalue is best, try first eigenvalue
             */
//            int n = ak_1.getCols();
//            Element ck = ak_1.getValue(n-1, n-1);
//            //Element ck = ak_1.getFactory().makeOne();
//            Matrix ckI = i.mult( ck );            
//            ak_1 = ak_1.sub( ckI );
            
            QRFactor fk = f.makeQRFactor(ak_1);
            
            
            qk = fk.getQ();
            rk = fk.getR();
            sk = sk.mult(qk);
            ak_1 = rk.mult(qk);
            
//            ak_1 = ak_1.add(ckI);
            
            /*
             * If we have converged bail early, note that I believe that
             * this test is too exact, the ones on the diagonal likely
             * converge faster than the off diagonals (zeros).  Also
             * there seem to be cases where the diagonal converges to 
             * a -1 instead of +1.
             */
            if( qk.isIdentity(e) )
            //if( qk.isDiagIdentity(e) )
            {
                //System.out.println("Q has converged to identity after "+k+
                //        " iterations");
                break;
            }
        }
        this.iterations = k;
        //System.out.println("Qk after "+ k+ " iterations converged "+this.converged()+ " = " +qk);
        // Last (previous) S has the orthonormal eigen vectors
        this.s = sk;
        this.q = qk;
        
        // The approximation of the eigen values should be on diagonal
        // The matrix should have small values below diagonal, e.g. U
        this.d = ak_1;
        
//        // Now determine the eigen vectors using LU to solve (A-eI)=0
//        Matrix identity = a.identity();
//        Vector zero = a.getFactory().makeVector( a.getRows() );
//        List<Vector> eigenVectors = new ArrayList<Vector>();
//        for (int k = 0; k < ak.getCols(); k++)
//        {
//            // Get the eigen value off the diagonal
//            Element ek  = ak.getValue(k, k);
//            Matrix a_ei = a.sub( identity.mult(ek) );
//            LUFactor f = new LUFactor(a_ei);
//            
//            Vector xk = f.solve(zero);
//            xk = xk.normalize();
//            eigenVectors.add(k, xk);
//        }
//        this.s = new Matrix(eigenVectors, a.getFactory());
//        
//        // Now no longer need ak, makeQRFactor diagonal
//        this.d = ak.diag();
    }
    
    public Matrix getA()
    {
        return this.a;
    }
    
    public Matrix getD()
    {
        return this.d;
    }
    
    public Matrix getS()
    {
        return this.s;
    }
    
    /**
     * Gets the matrix Q used in the iterations that should converge to 
     * the orthonormal, identity matrix I
     * @return 
     */
    public Matrix getQ()
    {
        return this.q;
    }
    
    public int getIterations()
    {
        return this.iterations;
    }
    
    public boolean converged()
    {
        return (this.iterations != this.maxIterations);
    }
}
